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We expand upon our previous analysis of numerical moving-puncture simulations of the 
Schwarzschild spacetime. We present a derivation of the family of analytic stationary l+log fo- 
liations of the Schwarzschild solution, and outline a transformation to isotropic-like coordinates. 
We discuss in detail the numerical evolution of standard Schwarzschild puncture data, and the new 
time-independent 1-l-log data. Finally, we demonstrate that the moving-puncture method can locate 
the appropriate stationary geometry in a robust manner when a numerical code alternates between 
two forms of l+log slicing during a simulation. 
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I. INTRODUCTION 

The binary black hole problem is a cornerstone prob- 
lem in gravitational theory: solving for the inspiral 
and merger of two black holes in full general relativ- 
ity connects the theory's strong-field regime with as- 
trophysics, and connects issues with the mathematical 
understanding of the theory with the emerging field 
of gravitational-wave astronomy. Solutions of the bi- 
nary black hole problem require numerical simulations, 
but stable, long-term numerical evolutions of orbiting 
black-hole binaries eluded researchers for four decades. 
However, after a number of insights and technical de- 
velopments two independent methods [U El E] were 
shown in 2005 to allow simulations of the last orbits, 
merger and ringdown of an equal-mass nonspinning bi- 
nary. One of these, the "moving puncture method" 
[21 13] , was quickly adopted by many research groups and 
has since been applied to many, more general, scenar- 
ios: unequal-mass binaries [H |51 151 17] , spinning binaries 

[S[3[ini[iii[ia[ni[ii[i3[ia[iii[iHi[iii[2ni'"2L 22, 23J, 

three-black-hole spacetimes [24l [25] and long simulations 
of (so far) up to ten orbits [ ^ ET l l ^ l^. 

At the technical level, the moving-puncture method 
consists of a seemingly small modification of the ear- 
lier "fixed puncture" method [30, 31 . However, in [32] 
(which we will refer to as Paper I), we found that the 
numerical slices behave in a way radically different from 
what was observed in previous fixed-puncture evolutions, 
and indeed in all previous numerical simulations of black- 
hole spacetimes. We developed a geometrical picture of 
the behavior of moving-puncture simulations, and our re- 
sults suggested a new paradigm for black-hole evolutions, 
centered around manifestly stationary representations of 
black-hole spacetimes and "asymptotically cylindrical" 
slices. This led to further investigations in [551 IMl 155]. 
In this paper we expand and discuss in more depth the 
results of Paper I, with particular reference to the asymp- 
totics of the initial and final states of a moving-puncture 



simulation 

The initial data in a typical moving-puncture simula- 
tion represent black holes using a wormhole topology: as 
we follow the coordinates toward one of the black holes, 
we do not reach the black hole's singularity but instead 
pass through a wormhole to another exterior space, and 
eventually find ourselves once again in an asymptotically 
fiat region. Data for N black holes can consist of -I- 1 
asymptotically fiat regions connected by A^ wormholes. 
In the puncture approach to constructing initial data 
[551 1571 1551 155] each "extra" asymptotically fiat region 
is compactified so that its spatial infinity is transformed 
to a single point ("puncture"), and the entire TV -I- 1- 
wormhole topology can be represented in a single three- 
dimensional space, M^. All of the black-hole singularities 
are conveniently avoided in this construction, and there 
is no need to "excise" any region when these data are 
used in a numerical simulation. 

This use of nontrivial topology to enforce the pres- 
ence of horizons can be regarded as a mere trick to con- 
veniently construct black-hole initial data. It comes at 
the expense of using the entire Kruskal extension to the 
Schwarzschild spacetime, part of which has no physical 
relevance in typical numerical evolutions. Our analysis 
suggests that a similar trick could be used to construct 
initial data that leave out most of the unphysical region; 
we will later refer to these as "trumpet" data. 

Standard puncture data are smooth over the entire 
space, except for a scalar function, the conformal fac- 
tor ^, which diverges as l/r near each puncture. One of 
the two innovations of the moving-puncture method was 
that it provides a method to stably evolve the conformal 
factor. In addition, the method uses gauge conditions 
(variants of "l+log" slicing for the lapse function [40] 
and P- freezing for the shift vector [SHUT]) that allow the 
punctures to move across the numerical grid. The result 
is that the punctures orbit each other and spiral inwards, 
as if the black holes were being represented by point par- 
ticles, and the plots of the "puncture tracks" shown in 
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many papers easily match our intuitive picture of ob- 
jects in orbit (see, for example, [21 |3j |26l |28] ) . Of course, 
the orbiting punctures are not point particles, they are 
the asymptotic infinities of wormholes, or at least they 
were in the initial data. Are these two extra copies of 
an asymptotically flat region of spacetime orbiting each 
other on the numerical grid? The answer turns out to be 
No, as we explain below. 

In Paper I we studied moving-puncture simulations of 
a Schwarzschild black hole, and found that the numer- 
ical points quickly leave the extra copy of the exterior 
space. Where the points near the puncture originally ap- 
proached a second asymptotically flat end, after a short 
time the grid points instead asymptote to an infinitely 
long "cylinder" with radius Rq = 1.31M. In addition, we 
solved the spherically symmetric Einstein equations for 
a stationary l-|-log-sliced spacetime (see also [42]), and 
found that the one solution with a lapse function that is 
non-negative everywhere agrees with that found by the 
niunerical code. In other words, there is one regular sta- 
tionary l-|-log solution, and the moving-puncture method 
quickly finds it. Whereas an embedding diagram of the 
initial data resembles a wormhole, we will refer to that 
of the late-time slices as a "trumpet" . This terminology 
becomes clear in Figures [l] and [2j 

This seemingly dramatic change in the appearance of 
the numerical slices is achieved by the gauge conditions. 
We start with momcnt-of-time-symmetry initial data for 
the Schwarzschild solution, choose the initial lapse a = 1 
and propagate using the l-|-log condition. So far the con- 
figuration is symmetric across the throat; we will call 
this "left-right" symmetry to be consistent with stan- 
dard Kruskal and Penrose diagrams, although we mean 
symmetry with respect to the upper and lower halves of 
Figure [T| So: the initial configuration is left-right sym- 
metric, and the slicing condition preserves this symme- 
try. This means that all the slices of the foliation will 
also be left-right symmetric and run from one spacelike 
infinity to the other. As such, there will not be a sta- 
tionary limiting slice; we will elaborate on this point in 
Section [IT] However, at late times a region of the slice 
to the right (and a corresponding region to the left) is 
asymptotically stationary. Now we have to consider the 
effect of the F-driver shift condition. If we use puncture 
data, this condition generates a very large shift, pointing 
to the right, near the puncture. This has the effect of 
dragging all the data points near the puncture into the 
region of the Schwarzschild solution between the hori- 
zons and finally onto the stationary part of the slicing. 
This happens extremely quickly. The closer the inner- 
most data point is to the puncture, the longer it takes, 
but it always happens. A description of this phenomenon 
was also given in [51] . 

There exists a true stationary l-|-log slice through the 
Schwarzschild spacetime. This is what we call the "trum- 
pet" . This is asymptotically flat at one end and cylindri- 
cal (of radius R ~ 1.31M) at the other. The asymptot- 
ically stationary part of the wormhole foliation asymp- 



totes to (part of) this trumpet. The closer the data points 
are to the puncture, the more of the trumpet we will fl- 
nally see. Further, we will see only the trumpet — there 
are no longer any grid points on the non-stationary part 
of the slice, or any of the slice on the left. We may say 
that the left half of the slice is grossly under-resolved 
(in fact, it is not resolved at all!), but the nature of the 
asymptotic l-|-log slice is such that a new asymptotics 
has formed, and the left half of the slice is causally dis- 
connected from the right: the non-resolution of part of 
the slice no longer concerns us. 

An example that closely mimics some of the behaviour 
described above was found in one of the first ever nu- 
merical simulations of a black-hole spacetime, as early 
as 1973, by Estabrook et al. [13]. They evolved the 
Schwarzschild solution with maximal slicing. The slices 
were reflection symmetric about a "throat" . At late times 
the left and right halves looked approximately station- 
ary. In fact the right half was translated to the right by 
the Killing vector while the left half was dragged to the 
left. The throat region approximated a cylinder of radius 
3Af/2, which grew linearly with time. 

Two phenomena associated with their simulations — 
"collapse of the lapse" and "slice stretching" — were re- 
curring topics in numerical studies of black-hole space- 
times for the next thirty years [HI US] HHJ H?]. How- 
ever, what we consider in this context to be the key fea- 
ture of their result, the formation of a cylindrical asymp- 
totics, received little, if any, attention before the work 
in Paper I. There we suggested that the existence of a 
time-independent representation of the foliation consis- 
tent with the new asymptotics, which a numerical code 
can find with the aid of an appropriate shift condition, 
is one of the keys to the success of the moving-puncture 
method. We also took the additional intuitive step of re- 
alizing that this time-independent representation could 
be utilized from the outset of a simulation, to allow the 
construction of fully time-independent trumpet puncture 
data. 

In this paper we extend the analysis that we provided 
in Paper I. In Section [IT] we discuss standard wormhole 
puncture data, which cannot be time-independent, and 
trumpet puncture data, which can. The canonical ex- 
amples are the maximally sliced solution [131 HH] and 
our stationary 1-l-log solution. Here we provide a simple 
derivation of the stationary l-|-log solution using a height 
function approach, and illustrate some of its features in 
a Penrose diagram. We also transform this solution to 
isotropic coordinates, providing puncture trumpet data 
that should be time independent in a moving-puncture 
simulation. 

In Section |III| we present our numerical method. We 
start with a brief description of the main features of the 
BSSN/moving-puncture system, and outline a procedure 
to construct Penrose diagrams from numerical simula- 
tions of the Schwarzschild spacetime. Scction |TV] contains 
the first set of numerical results: a detailed study of a 
moving-puncture evolution of wormhole puncture initial 
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data. We show how quickly the numerical data approach 
a trumpet geometry, how they behave during the transi- 
tion, and check the accuracy with which they approach 
the analytic stationary solution. It is important to em- 
phasize that although the numerical solution is stationary 
in the sense of coordinate-independent functions (for ex- 
ample, a plot of a vs i? or of Tt{K) vs a), the numerical 
coordinates may still exhibit some drift. 

This is emphasized in Section [V] where we evolve 
the trumpet puncture data produced in Section [ll] We 
demonstrate that these are indeed time independent, up 
to small numerical errors. Furthermore, if we alternate 
between two variants of H-log slicing during an evolu- 
tion, we show that the numerical data can alternate be- 
tween two stationary solutions, but that some coordinate 
drift develops in the process. This drift is minimized by 
choosing r; = in the F-driver shift condition. 

Finally, we conclude with some remarks about the po- 
tential usefulness of trumpet puncture initial data for 
black-hole binaries, which will be the subject of future 
work. 



II. WORMHOLES AND TRUMPETS: 
ANALYTIC TREATMENT 

We begin by summarizing our notation and approach, 
which is based on the standard "3-1-1" space-|-time 
splitting of Einstein's equations [W. We introduce 
the Schwarzschild metric in Schwarzschild coordinates, 
progress to isotropic coordinates, which are better 
adapted to the standard "puncture" method, and then 
on to a derivation of the "trumpet" solution previously 
introduced in Paper F Our focus here is on solutions that 
are time-independent, because even in the black-hole bi- 
nary problem we want to choose coordinates that min- 
imize the gauge dynamics; the only dynamics we really 
want to see in a numerical code are physical dynamics. 



sign convention (SDJIHT], not that given in Wald It 
proves convenient to also split the extrinsic curvature into 
its trace K and tracefree part Aij, i.e., 

K,, = A, + ^7,,if. (3) 

We further decompose these data with respect to a 
conformal metric, 7^^ , as follows: 

7y = (4) 

Aij — t(j ^Aij, (5) 

K ^ K, (6) 

= P'- (7) 

Note that in the standard conformal decomposition of- 
ten used to solve the constraint equations [53] , one 
chooses p = —2, but for the BSSN decomposition used 
in the moving-puncture method, p ~ A. Note also that 
the trace of the extrinsic curvature and the contravari- 
ant components of the shift vector are unchanged by the 
conformal rescaling. The lapse function can also be trans- 
formed, but this is not necessary for the present study. 
We could also have made other choices of the conformal 
weights, but we are not interested in those here. The 
complete data on one time slice are {-jijjKij), or equiv- 
alently {i>,jij, Aij, K). 

For the moment we will not worry about the details 
of the constraint or evolution equations. In this section 
we will be interested in finding time-independent repre- 
sentations of the Schwarzschild solution, for which the 
constraint equations are already satisfied, and the time 
evolution of the data is trivial (i.e., they do not change in 
time). Although we will make comments about numeri- 
cal simulations in this section, we will postpone concrete 
details until Section |III[ when we study numerically the 
time development of our data. 



A. 3+1 decomposition and variables 

We start by making a space plus time (3-1-1) decompo- 
sition of the spacetime metric, 

ds^ = -a^dt^ + -i^j{dx' + (3'dt){dx^ + (3^dt). (1) 

The three-dimensional metric of a t = const slice is de- 
noted by . The lapse function a gives the proper time 
between the slice at time t, and the next slice at time 
t + dt. The shift vector (3^ prescribes how the coordinates 
shift between the two slices. The complete data on one 
timeslice are given by 7^ and the extrinsic curvature 

= 1^ ( + V, A - da,j ) , (2) 



B. Wormhole puncture data for a Schwarzschild 
black hole 

The Schwarzschild metric in Schwarzschild coordinates 

is 

ds^ = - fdT^ + f-^dR^ + R^dn^, (8) 

where / = 1 — 2M/R. Throughout this paper R 
and T denote the Schwarzschild radial coordinate and 
Schwarzschild time. The surface R = 2M is the event 
horizon, i? = is a physical singularity, and i? ^ 00 is 
spatial infinity. 

If we make the coordinate transformation R — tp'^r, 
where 



where denotes covariant differentiation with respect to 
the spatial metric jij, and we have used the MTW/ADM 



(9) 



FIG. 1: Embedding diagram of a two-dimensional slice (T = 
const., ^ = IX 1 2) of the extended Schwarzschild solution. The 
distance to the rotation axis is R. A wormhole with a throat 
aX R = 2M connects two asymptotically flat ends. 



the Schwarzschild metric becomes 

= dT^ + i>\dr^ + r^d^''). (10) 

These are called (quasi or spatially) isotropic coordinates. 



Topologically, the constant-T slices are 



X 



\ 



{0}, while geometrically (measuring proper areas or the 
Schwarzschild radius R) the slices are wormholes. The 
isotropic r does not reach the physical singularity at i? = 
0. For large r we see that R oo, but for small r we 
see that once again i? — > oo. There is a minimum of 
R = 2M at r = M/2. We now have two copies of the 
space outside the event horizon, R > 2M, and the two 
spaces are connected by a wormhole with a throat at 
R — 2M (Figure [ij . This wormhole picture of a black 
hole forms the basis of the initial data used in current 
moving-puncture black-hole simulations. The point r = 
0, which represents the second asymptotically flat end, is 
referred to as the puncture. 

In terms of the conformal 3-1-1 quantities introduced 
earlier, the Schwarzschild metric in isotropic coordinates 
is 



= 1- 

j = , 

K ^ 0. 

The lapse and shift are 



M 
2r' 



1 - 

IT 

0. 



M 
2r 
M ' 
2r 



(11) 

(12) 

(13) 
(14) 



(15) 
(16) 



If we choose ( 11 1 - ( 14 ) as our initial data, and propagate 



the data with the lapse ( 15 ) and shift ( 16 1, then the data 



FIG. 2: Embedding diagram of a two-dimensional slice (t = 
const., ^ = n/2) of the maximal solution (171 - (201. The 
distance to the rotation axis is R. In contrast to Figure [l] 
there is only one asymptotically flat end. The other end is an 
infinitely long cylinder with radius Rq = 3M/2. 



As trivial as the time development of these data is, 
it is difficult to reproduce numerically in a standard 3D 
black-hole evolution code. The reason is that most codes 
are not stable when the lapse is negative, which it is here 
for r < M/2. In a numerical code we prefer to use a lapse 
that is always positive, or at least non-negative. Unfortu- 
nately, it is not possible to construct a time-independent, 
maximal slice of Schwarzschild with two asymptotically 
flat ends and an everywhere non-negative lapse, no mat- 
ter what shift conditions we employ [53J [55J [35] . Put an- 
other way, a maximal- or 1-1- log-slicing evolution with two 
asymptotically flat ends and with a non-negative lapse 
cannot reach a stationary state. We will now see that 
the way around these earlier results is to give up one of 
the asymptotically flat ends. 



C. Trumpet solution — maximal slicing 

Since maximal initial data with two asymptotically flat 
ends cannot be time independent, we seek an alternative. 
One alternative is "trumpet" data. Consider the maxi- 
mal slice of the Schwarzschild spacetime |13J H5] 



- 1- 



2M 

diag(2C/i?^ -C/R^,-C/R^), 
aC 



2M C2 



1 



R R^' 



(17) 
(18) 
(19) 

(20) 



will remain unchanged: they are time-independent data 



with C = 3y3MV4 and R € [1.5M,oo). We see that 
a > for this domain, and at Rq = 3M/2 the lapse goes 
to zero, as does /3^', while the spatial metric diverges as 
R Rq, and so the proper distance from Rq to any 
R > Rq is inflnite. As we approach Rq the timeslice 
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becomes an infinitely long cylinder of radius 3M/2. This 
behaviour is illustrated in Figure [2] and we will refer to 



these data as "trumpets" . These data, like pT 
with (15) and (16 1, are time independent 



(14) 



With the benefit of hindsight we see that, rather than 
adopting the wormholc puncture data ( 11 1 - (14) for nu- 



merical evolutions, we might be better off transforming 
(17) - (20) to isotropic-like coordinates (such that r — 
corresponds to i? = 3M/2), and using those data instead. 
Such a transformation was calculated numerically in |33j . 
and it was shown that these data can indeed be evolved 
stably, and are time-independent up to (convergent) nu- 
merical errors. An analytic transformation to isotropic 
coordinates was given in [57] . 



D. Trumpet solution — 1+log slicing 



The data ( 17 1 - ( 20 1 are maximally sliced, K = 0. In a 



numerical simulation, we must solve an elliptic equation 
at each time step to find the appropriate lapse function 
that maintains maximal slicing (for these data we could 
assume that the lapse remains constant, but this will not 
be true in general black-hole simulations). This is com- 
putationally expensive (i.e., slow), and it is therefore cur- 
rently more practical to choose a slicing condition so that 
the lapse can be calculated from an evolution equation 
in the same way as all of the other dynamical variables. 
One such slicing condition is 



dta — ~naK, 



(21) 



where n is some constant, usually chosen to be ti = 2. 
We see that \i K = 0, then the lapse does not evolve, 
and this condition will maintain maximal slicing for our 
trumpet data. This slicing condition is part of a class 
of conditions called "1-t-log" slicing [3Cr. The condition 
(21 1 is actually not a "geometric" slicing condition in the 



sense that the slicing resulting from this condition also 
depends on the shift /3% i.e., on the spatial coordinates 
[40, 58J. Furthermore, in binary simulations it has been 
found that unphysical gauge modes arise when (21) is 
used [59], and a preferred choice is 



{dt - (3''di)a = -naK. 
This is equivalent to 

Cf,a — ~nK, 



(22) 



(23) 



where n is the unit timelike normal to the slice and C is 
the Lie derivative. For the rest of this article, whenever 
we refer to l-l-log slicing, we will mean (22). We will 



refer to Eqn. (21) as "asymptotically maximal slicing" 



because, if it leads to a time-independent geometry, then 
that geometry will be maximally sliced. 

Since the slicing condition ( 22 ) has proven rather ben- 



eficial in black-hole simulations, we would like to know 
what the corresponding time-independent Schwarzschild 



data are. A stationary l-|-log solution was found in Pa- 
per I. An earlier result on stationary l-|-log slices can be 
found in [12], however its exact relation to the present 
discussion of punctures remains to be understood. Here 
we will present one (of many possible) derivations of the 
solution of Paper I. 



1. Height function derivation 

We begin with the Schwarzschild metric in 
Schwarzschild coordinates, Eqn. ([s]), and introduce 
a new time variable, t, related to Schwarzschild time by 
a spherically symmetric height function, h{R), as in [60] . 



t -h{R). 



(24) 



The lapse, shift and future-pointing normal vector in the 
new coordinates are 



/ 



1 - /2/l'2 ' 



Ph'^ - 1 ' 



= (-a, 0,0,0), 



(25) 

(26) 
(27) 



Here again / = 1 — 2M/R, and we have introduced the 
notation h' = dh{R)/dR. We note that only the deriva- 
tive of the height function appears in the new metric. 

It is also useful to note that for any stationary spheri- 
cally symmetric spacetime the extrinsic curvature is given 

by 



/3' 



Kee = i?/3, 



(28) 

(29) 
(30) 



where (3 = \J PiP^ and /?' = d(3/dR, and the trace of the 
extrinsic curvature K = Kl is given by 



(31) 



For the Schwarzschild solution we also have the relation 

2M 



a' -P' = 1- 



R 



(32) 



It will simplify the following calculations if we write h' 
in terms of the lapse as 



h' 



fa ' 



and the shift is now written as 



/3« = a^a^ ~ f, 



(33) 



(34) 



where we have chosen the sign of the root to give a non- 
negative shift. 
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Following [BU], the trace of the extrinsic curvature, K, 
can be written as 



K 




dniR'h'fa). 



(35) 



In turn, our 1+log gauge condition can be written for a 
time-independent solution as 



(3 d^a — naK. 



(36) 



Combining Equations ( 35 1 and ( 36 ) and solving for a' 



we arrive at the following first-order differential equation 

for the lapse. 



a 



n{3M - 2R + 2Ra^) 
' R{R - 2M + nRa - Ra'^) ' 



(37) 



The solution of this equation is an implicit equation for 
the lapse, 

- [3 In i? + ln(2Af - i? + Ra^)] - a = constant, (38) 



or, more conveniently. 



= 1 - 



2M 



C(n)2e2"/" 



(39) 



We will now determine C(n). The equation (37 1 is 
singular when the denominator is zero. We search for a 
C(n) such that the numerator and denominator are zero 
at the same point, and the equation remains regular (the 
singularity in the equation signifies a transition between 
an elliptic and a hyperbolic problem as discussed in |35p. 
The numerator is zero when 



a = ±a/1 



3M 
2^' 



(40) 



We want solutions with an everywhere non-negative 
lapse, so we choose the positive root. Substituting this 
into the denominator of ( 37 ) we find that the denomina- 



tor is zero for a particular value of R, 



Rr — 



All? 



(41) 



where we have chosen the positive root to ensure that Rc 
is always positive. This allows us to evaluate the lapse 
at that point as 



V4-f 9r 



3n 



The value of the constant C{n) is now given as 



(42) 



C-in) = (5!i±^i+M)!e--./". (43) 
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FIG. 3: Location ol the throat, _Ro, as a lunction of the coef- 
ficient n in the slicing equation (36 1. In the limit n — > oo the 



stationary slice is maximal and the throat is at Ro = 1.5M, 
indicated in the figure by a horizontal line. For small values 
of n, the throat approaches the singularity. 



We now have the 1-1- log trumpet solution, given by 



(39 1 with (42) and (43 1 for the lapse, from which we 



calculate /?^^ from ( 34 1 and 7^7? — , and the extrinsic 



curvature is given by ( 28 1 - ( 30 1 



Three special cases of the constant n deserve immedi- 
ate comment. When n = 2 we have the form of l-|-log 
slicing commonly used in numerical simulations of black- 
hole spacetimes. The results presented in [351 [Ml [3S] 
pertain to this case, and we have Rc — 1.54M, C^(n) = 
1.55431M'', and the location of the throat (the small- 
est real root of ^ with a = 0) is i?o = 1.31241M. 
The horizon is located at a{R = 2M) = 0.376179, which 
differs by over 20% from the common "a = 0.3" rule-of- 
thumb estimate of the horizon location. 

In the limit n —t Q, Eqn. (36) is degenerate, and two 
solutions exist. In one, — 0, and we recover the stan- 
dard Schwarzschild metric in Schwarzschild coordinates. 
In the other, dj^a = 0, and so a = 1 everywhere. As we 
approach n — > 0, the radius of the throat i?o approaches 
the singularity as i?o ~ e-2/(3"V(2^/^«)- 

The third case of interest is the limit n ^ 00. Now 
the only physically meaningful solution is that if = 
everywhere, i.e., maximal slicing. We have C(oo) = 
3\/3Af^/4, and the radius Rc and the throat coincide 
at Rc = Ro — 3M/2. This is the cylindrical maximal 



slice (171 ~ (20) 



2. Penrose diagrams 

One advantage of working in spherical symmetry is 
that redundant coordinates may be suppressed and we 
can visualize the way the spacetime is sliced on two- 
dimensional diagrams, such as the Carter-Penrose dia- 
gram. 



In order to do this for a given solution (39), we first 



integrate the height function h using (33) to obtain T 



for every R. The singularity at i? = 2M can be handled 
(at least numerically) by introducing a different quan- 
tity, such as e~'', around the horizon. The undetermined 
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FIG. 4: The Penrose diagram of the shoes defined by the 
stationary solution of the 1+log condition with n = 2. Every 
slice approaches along the curve R — Ro ~ 1.31Af and 
spatial infinity i'^ along a curve of constant T. The slices are 
displayed in time steps of 4M. 



value of t in Eqn. (24), which can be interpreted as the 
constant of integration, expresses the fact that we do not 
calculate a single slice but a foliation of the Schwarzschild 
spacetime. As expected, the slices are related to each 
other by sliding along the Killing vector field dx- 

From the coordinates {R, T) along one slice we trans- 
form to Kruskal coordinates (u, v) by either 



for i? > 2M or 



R 




2M 


- 1 


R 




2M 


- 1 


1 


R 




2M 


1 


R 




2M 



1 e«/(4M) ^-^^ 



3^/(4^^) sinh 



T 

4M' 
T 

4M' 



T 
4M' 
T 

4M' 



for R < 2M. Compactifying the result via 
u±v = tan(C/±y), 



(44) 
(45) 

(46) 
(47) 

(48) 



where U and V are the abscissa and the ordinate of 
the Penrose diagram, yields the picture displayed in Fig- 
ure [4] Note that we chose non-negative shift which 
also determines the sign of h (as well as u and v in Eqns. 
( [44| and (45l) . The opposite choice would lead to slices 
that are mirror images of those in Figure [4] connecting 



to i 



3. Trumpet data in isotropic coordinates 

We have now derived the stationary l-|-log solution. 
In Paper I we compared this solution with the late-time 
data from moving-puncture evolutions of wormhole punc- 



would also like to put this solution into isotropic-like co- 
ordinates, as was done for the maximal trumpet data in 
[551 [57] . This provides a good test-case for numerical evo- 
lution codes, and could be a starting point for generating 
trumpet data for black-hole binaries. 



The implicit nature of the solution ( 39 ) makes it dif- 
ficult to analytically construct the transformation to 
isotropic coordinates. However, solving (39 1 for a func- 



tion R{a) leads to four roots of a fourth order polynomial, 
which, if chosen appropriately, represent the analytical 
solution. In this section, R should always be understood 
as this function of a, whereas a becomes the indepen- 
dent variable that parametrizes the spatial dependency. 
Apart from that our approach is similar to the one used 
in [ST . We note that the ^rr component of the station- 
ary l-|-log metric can be related to the 7rr in isotropic 
coordinates by 



IRR 



dr 
dR 

dr 
dR 



We therefore find that, using R — tp'^r and jrr — a 



dr 
OR 



r 
oR' 



and a relation between the isotropic coordinate r 
Schwarzschild R may be found by either 



exp 



exp 



1 dR 
otR da 

1 Ini? 



{a) da 
da 



(49) 
(50) 

(51) 
and 

(52) 
(53) 



where the last equation is obtained by integration by 
parts and the upper integration limit is chosen such that 
r— >_RasQ;— >l,i.e., towards spatial infinity. Both in- 



tegrals ( 52 1 and ( 53 1 diverge as a 

;ly, and can be integrated numerically to arbi- 



0, but (52 1 diverges 

less strong! 

trarily small a with sufficient accuracy to produce data 
suitable for a numerical evolution. On the other hand, 
( 53 1 has the attractive property that as a ^ 1 the factor 
RP"" gives the asymptotic behaviour that we wish. In 
practice, we choose a point a^ = 0.1, and for a < a^ use 



ture data, which we will discuss again in Section III We 



r(a) = i?(a,)(i/"=)exp 
where 

Co = 
For a > as we use 

r(a) = i?(a)(i/")exp 



1 dR 

Q aR{a) da 

^ InR(a) , 
— da. 



(a) da — Co , 
(54) 

(55) 



In R{a) 



da — C'o 



(56) 
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Once we have the solutions R{a) and r{a), we may con- 
struct r{R) to whatever accuracy is desired, and then 
transform our data to isotropic coordinates via 




K 



RR- 



(57) 
(58) 
(59) 



Note that the singularity in the conformal factor is 
now milder than in the standard puncture case, where 
ip ~ 1/r, while here 



il){r) 




(60) 



This fact is intuitively expressed in the embedding dia- 
grams in Figures [l]and|2] the expansion of the wormhole 
geometry as compared to the trumpet geometry leads 
to a stronger singularity of the conformal factor. Since 
wormholes and trumpets both allow a representation on 
K'^ where a coordinate singularity at r = is absorbed in 
the conformal factor ip, we refer to both cases as punc- 
tures. 

All of the numerical calculations described here were 
performed with Mathematica, and the data output as ta- 
bles of physical quantities parametrized by the isotropic 
coordinate r. The data files were then read into our 
full 3D code [HU [55], where they were transformed to 
Cartesian coordinates and interpolated onto the numeri- 
cal grid. Examples of the resulting data for a, (3^ and K 
are shown in Figure [5] The time independence of these 
data will be explicitly demonstrated in numerical evolu- 
tions in Section fVl 




III. NUMERICAL SIMULATIONS 

We now turn to the numerical evolution of wormhole 
and trumpet puncture data for the Schwarzschild space- 
time. We start with a brief description of the moving- 
puncture method and our numerical techniques, as well 
as a summary of our procedure for analyzing our results 
(including the construction of Penrose diagrams), before 
finally presenting our numerical results. 



A. The BSSN/moving-puncture system 



The 3-1-1 decomposition provides evolution equations 
for the spatial metric ^ij and extrinsic curvature Kij |49j . 



The BSSN reformulation consists of rewriting the evolu- 
tion equations in terms of conformally rescaledvariables 

-Aijy where we now use p 
additional variable is introduced: F* : 



4 in ([5]), and an 



FIG. 5: The lapse a, a;-component of the shift vector /3^, and 
K for time-independent 1-l-log data in isotropic coordinates. 
The data are shown along the a::-axis. 



When we deal with puncture data, the conformal fac- 
tor ip diverges at each puncture, and this is handled in 
the moving-puncture modification [21 [3] of the BSSN sys- 
tem by replacing ip with either the variable (p = Inip or 
X = V'"'' (or X = V'~^ [SSj), and one of these quanti- 
ties is evolved instead. The details of the BSSN system 
are given in [MJ We use the BAM code, and pro- 
vide details of our implementation of the BSSN/moving- 
puncture system in jSTj . 

Given evolution equations for the variables 
{jij, Aij, K^T"^} and (p or x, and some initial data, 
we also need to choose a lapse and shift during the 
evolution. We have already discussed the H-log evolu- 
tion equations (21 1 and (22 1 for the lapse function; the 



stationa rv sli ces for these slicing conditions are given in 
Section IIP For the shift vector we use the F-driver 
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condition [511111] , 



J u dr 



dtB' 



(61) 
(62) 



This shift condition is crucial to the behaviour of the 
moving-puncture system. In black-hole binary simula- 
tions it generates a shift that moves the punctures around 
the grid on trajectories that match very well the motion 
that would be seen from infinity [27l|32]. This shift also 
allows the wormhole puncture data to evolve to the sta- 
tionary l-|-log geometry, as shown in Paper I, and as we 



will see again in Section IV As we said in the Introduc- 
tion, the approach to the "puncture geometry" requires 
a dramatic stretching of the coordinate representation of 
the slices. We will illustrate this extreme behaviour in 
detail in the coming sections, but it can be seen most 
immediately by calculating the norm of the shift vector, 
= 7jj/3'/3-', during the first few M of evolution: al- 
though /?' remains finite, and in fact goes to zero as we 
approach the puncture, diverges. This can also be 
seen analytically [66]. The wormhole slice is stretched 
such that all of the numerical points extremely quickly 
leave the part of the slice that cannot be stationary, and 
the part of the slice in the second copy of the exterior 
space, and the points relax onto the stationary 1-1- log 
slice. This could not happen with a zero shift, and of 
course would not happen with an arbitrary shift. There 
may be a large class of shift conditions that produce the 
same effect, but the first to be found, and the one that is 
standard in moving-puncture simulations, is the F-driver 
condition (611 - (62), and variants. 



B. Penrose diagrams from numerical data 

A convenient way to view the numerical evolution of 
our data is to represent them in a Penrose diagram. 
Given the Schwarzschild coordinates R and T of a given 
point, it is straightforward to calculate the corresponding 
point on a Penrose diagram. From our numerical data we 
can easily calculate the Schwarzschild coordinate, R{r, t); 
see Section |IV| and Eqn. [73] Each point on the initial 
slice is at a constant Schwarzschild time T(r, 0) = Tq, 
which we are free to choose, but T(r,t) is not known. 
Although a variety of ways to compute an appropriate 
coordinate representation suggest themselves, it actually 
turned out to be somewhat tricky to find one that works 
reliably with our numerical data, see also the discussion 
in [67] where the transformation to Kruskal coordinates 
is implemented. Care has to be taken near the horizon 
at i? = 2M, and in our case some issues arose at mesh 
refinement boundaries. 

An outline of the procedure that we settled on is as 
follows. The method is also illustrated in Figure [6] We 
first choose a numerical point tq far from the black hole 
(in practice a numerical point that does not pass through 




FIG. 6; A sketch of the method we use to obtain T{r, t). The 
lines are drawn from actual data with ro = 3.25M for the 
integration in time and t — lOM for the integration in space. 
Only a subset of the grid points is displayed. 



the horizon during the evolution). As we have already 
said, on the initial slice that point has Schwarzschild time 
r(ro,0) = Tq. We use a differential equation for T — 
dT/dt to integrate forward in time through our numerical 
data and produce T(rQ,t). Then, on each time slice, t — 
ti, we use a second differential equation, this time for 
T' — dT/dr, to integrate across the slice and calculate 
T{r,ti). The equation we use is badly behaved near the 
horizon, and for a set of five points across the horizon we 
integrate instead a differential equation for u' — du/dr, 
where u is the Kruskal coordinate defined by (44) and 



(46) 



We derive the required differential equations by consid- 
ering the transformation of the spacetime metric between 
Schwarzschild and numerical coordinates. A differential 
equation for T can be found using the transformation 



gtt 



We therefore find 



dT 
~dt 



9tt 



dR 
~dt 



9RR- 



(63) 



2M\ 
If) 



gtt 



(64) 



This expression is valid for R > 2M and m > 0, and 
we choose rg — 3.25M, where this is always true. The 
metric component gtt is given by gtt — — (a^ — (3iP^). 
Using ( 64 1 we can integrate the Schwarzschild time along 
the constant — 3.25M curve for the duration of the 
simulation. For the numerical coordinate tq = 3.25M we 
now know R and T throughout the numerical evolution. 

A differential equation for T' on a slice of constant 
numerical coordinate time can be found by transforming 
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the spatial part of the metric, 



dT 
dr 



dRX 
dr J 



9rr 



^ rj,l2 ^ 7rr - R' 9B.R 

9rr 



R'2 h 2M\ 



(65) 

(66) 
(67) 



Note that R and R' can be readily computed from the 
numerical data, since we are dealing with data on one 
numerical timeslice. In order to integrate (67 1, we must 



take a square root and choose the sign such that the slices 
go smoothly through the horizon, and to define which 
ends of the computational domain belong to which side 
of the Penrose diagram. The result is 



r 




2Af\ 



(68) 



Equation ( 68 1 is integrated along the entire numerical 



slice, except at the points near the horizon, where T' 
is singular. We overcome this difficulty by instead inte- 
grating the Kruskal coordinate u through five points that 
cross the horizon. A differential equation for u' can be 
found by transforming the {r,t) coordinates to {u,R): 



du 
dr 



9u'b 



du OR 

^~^~a~9uR 

or or 



dR^ 



9RR- 



(69) 



The metric components guu, guR and guB, can be writ- 
ten in terms of only u and R by making use of the 
Schwarzschild metric in Kruskal coordinates, g^^: 



V 

guu 

9uR 

9rr 



dv dv 



\2M 

dv 
du 



1 e-R/(2^-f) (^; > 0), (70) 



gv 



1 



32M 
~~R 



43 



K 



dudR^"" M' 



dv \ 

dR J -^"^ 



K 



4m 1 
M 
R 1 



R/(2M) 



,-R/{2M) 



(71) 



Producing an equation for u' from ( 69 1 once again re- 
quires the appropriate choice of sign for a square root. 
The final result is 

= — (-9urR' + J {9uRR'f + 9uu {irr " R'^9Rr)] 
9uu \ / 

(72) 

Once equations ( 68 1 and ( 72 1 have been integrated along 



a slice, a constant of integration is chosen to give the 
value already calculated at rg = 3.25M using (64 1, and 
a choice of Tq for the time on the initial slice. 



IV. THE NUMERICAL TRANSITION FROM A 
WORMHOLE TO A TRUMPET 

A. Setup 

Simulations were performed with computational grids 
consisting of eight nested cubical boxes. The boxes are 
labeled by Z = 0, 1, . . . , 7. Each box contains points, 
and the grid-spacing in each box is hi = H/2'-, where 
H is the grid-spacing of the largest box. We denote 
hmin = ^;,„ai ■ The simulations were performed in only 
one octant of a three-dimensional Cartesian grid (the 
data in the other octants being easily deduced by exploit- 
ing the known spherical symmetry of the Schwarzschild 
spacetime). The grid points are staggered across the 
coordinate axes, such that the boundaries of the cubi- 
cal domain are the six planes defined by a; = /imm/2, 
y = z = hrmnf^, x = X, y = Y and z ^ Z, 

with X = Y = Z = 192M, which define the stan- 
dard "outer boundary" of the computational domain. 
Three simulations were performed to make a convergence 
series, with N = 64,96,128, and H = 6M,4M,3M, 
h„nn = Af/21.33,M/32,M/42.67. Note that these grid 
configurations are the same as used for the single black- 
hole tests in [61j . 

The numerical simulations discussed in this section 
start with the initial data (111 - (14 1, and initial shift 

= and initial lapse a = 1. A number of additional 
simulations were performed with initial lapse a = ijj^"^ 
for comparison. 

The results of the numerical evolution are shown in 
Figures [7] to [18] The data for these figures were pro- 
duced by interpolating numerical data onto the x co- 
ordinate axis. The spherical symmetry of the solution 
allows us to analyze the simulation using these data 
only. The Schwarzschild radial coordinate can be cal- 
culated by relating the numerical spatial metric to the 
Schwarzschild metric ([sj. We know that — R^, 
and because of the spherical symmetry we also have 
-fee = {dx'/de){dx^/de)-f,j = {dx' /d9){dxj /de)ip%j. 
Along the x-axis, 9 = 7r/2, and so dx/dO = dy/dO = 
and dzjdd = x, and we have 



R |l; = z=0 



(73) 



B. Early-time behaviour 



Figure |7] illustrates the main global feature of the time 
development of the numerical slices, namely the transi- 
tion from wormhole to trumpet asymptotics. The upper 
figure shows the proper distance of a given point from 
the horizon at i? = 2M versus that point's Schwarzschild 
coordinate R. The thick line indicates the t = data. 
Initially R — 2M is the minimal surface, and the data 
contain two copies of the space outside the horizon. This 
is the initial wormhole that we have referred to several 
times. If we consider a surface of revolution around the 
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FIG. 7: The top two figures show the proper separation from 
the (outer) horizon versus the Schwarzschild coordinate R. 
The left panel shows the slices a,t t = 0, l,2,3Af, and the 
right panel shows the slice at t = 50M. The final numerical 
slice terminates at ii ~ 1.31M. The vertical line indicates 
the horizon a,t R = 2M , and the six points represent x/M = 
1/40,1/20,1/8,2,5,8 on each shce. The lower two figures 
show Schwarzschild R versus numerical r for the same points 
at the same times. The horizontal lines show R — 2M and 
R = 1.31M. 



proper time axis, we obtain just another version of the 
standard picture of a spacetime wormhole, as shown ear- 
lier in Figure [T] We show the data only up to coordinate 
R = lOAf ; the upper and lower Hues in the figure quickly 
asymptote to ±45 degrees, and do not add extra informa- 
tion to the figure. If we follow the upper line outwards, 
we move further from the origin (the puncture) , and the 
plot ends at r = 8.25M. If we follow the lower line out- 
wards, we move closer to the puncture, and in this case 
the plot ends at r « M/40. The corresponding picture of 
R in terms of coordinate r is shown in the lower panels 
of Figure [Tj 

At early times two notable things happen. First, the 
minimal surface shifts to i? < 2M, and the numerical 
domain contains two surfaces with R — 2M, which we 
will call the "inner" and "outer" horizons. The proper 
distance in Figure [7] is with respect to the outer hori- 
zon. Second, the points in the lower right part of the 
figure rapidly move to the left. In other words, the 
Schwarzschild R corresponding to those points rapidly 
decreases. The numerical points don't "go" anywhere, of 
course; they are points on a fixed grid. But their location 
in the Schwarzschild spacetime does change, and quickly. 
Within only 3A/ all of the points have passed the inner 
horizon, and it has ceased to be part of the numerical 
domain. 



At later times {t > 40M) the points close to the punc- 
ture settle at a constant value of R, which we will soon 
see is close to Rq = 1.31 A/, the location of the throat 
in the stationary H-log solution derived in Section [IIP [ 
This is shown in the upper right figure. At first sight this 
does not correspond to the cylindrical asymptotics shown 
in Figure [2] the "cylinder" looks too short. The reason 
is that the spatial metric now diverges more slowly as 
we approach the asymptotic region, and so a point ini- 
tially "close" to the second asymptotically flat end was a 
larger proper distance from the outer horizon than it is 
now that it's "close" to the cylinder. 

The lower two panels in Figure [7] give a complementary 
picture. These plot Schwarzschild R versus the numerical 
coordinate r to illustrate directly how the Schwarzschild 
R of each grid point changes. The behaviour at early 
times in terms of R{r, t) is also shown in Figure [s] which 
shows several contour plots of R (indicated by variations 
in colour) at different values of coordinate r and time t. 
These figures illustrate many of the main features of the 
early-time evolution of the wormhole puncture data. We 
can clearly see that the horizon, initially at r = M/2, 
splits into two copies, and also that values of i? < 2M, 
which are not present in the initial data, appear as time 
progresses. We can also see that for larger r the value of 
R changes very little, although it does decrease for all r. 

The behaviour of grid points close to the puncture is 
shown in Figure [9j This plot uses data from a simulation 
that used extremely high resolution at the puncture: 16 
nested boxes, each containing 64"^ points, with a coarsest 
resolution oi H = AM and a finest resolution of h,nin = 
Af/8192; the simulation was run for t — 5M, in order to 
obtain the results displayed in Figure [9] 

As shown in Paper I, the numerical data become dis- 
continuous across the puncture. This means that finite- 
difference derivatives (which are used to calculate many 
quantities in the BSSN evolution) will have even worse 
discontinuities, and the numerical method cannot con- 
verge for points that are within a stencil-width of the 
puncture, which for the simulations discussed here in- 
cludes the two grid points closest to the puncture. Fortu- 
itously the nature of the BSSN/moving-puncture system 
is such that these errors do not seem to propagate out 
from the puncture, and clean convergence can be seen 
up to the last few grid points. (This is shown in Fig- 
ures 2 and 3 of fSI].) In the upper panel of Figure [9] we 
show the time development of R for the third, fourth, 
and fifth closest points to the puncture of our extremely 
high-resolution simulation, at r = {3, 4, 5}M/8192. Al- 
though the points closer to the puncture are not expected 
to show clean convergence, we see in the lower panel of 
Figure [9] that they display similar behaviour. 

Initially the points in the upper panel are at R— (1 + 
2M/rfr = {684, 513, 410}M. The value does not change 
significantly for the first M/2 of the simulation, but then 
quickly deceases. 

How quickly are our numerical slices flung out of the 
second copy of the exterior space? Figure 10 shows 
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FIG. 8: Contour plots showing the behaviour of R(r, t) on the numerical grid for the first 20M of a simulation. The axes of 
the figures give r and t in units of M for three different coordinate ranges. The contour-labels and the colours indicate the 
value of R. We clearly see that R changes rapidly for a given value of the r-coordinate for the first few M of evolution before 
settling down to approximately stationary values at later time. Points with r > M/2 accelerate towards the black hole {R 
decreases) before settling down again at some smaller value of R. As the slice moves towards the Schwarzschild singularity, the 
horizon, initially at r = M/2 and R — 2M, splits into two copies, and values of i? < 2M, which are not present in the initial 
data, appear as time progresses. The initial slice covers the interior of the black hole from R = 2M to 71 = +00 for r = M/2 
to r = 0. This region is quickly squeezed towards zero extent in r, as the lines of constant R in the lower left of the panels 
indicate. Note that for given r, the coordinate motion R{r, t) is not monotonic in t, but R{r, t) approaches its asymptotic value 
via a damped oscillation, see also Figure [Tl] 



600 
500 







\ \ 


— r = 3M/8192 




r = 4M/8192 




r=5M/8192 















1 2 3 4 5 

Time (M) 



2000 



S. 1000 

a! 









r = M/S]92 
r= 2M/S]92 


N \ 


r = 3M/8192 


^ S \ 













12 3 4 

Time (M) 



FIG. 9: Evolution of R for fixed r for the first 5M of an ex- 
tremely high-resolution simulation. Top panel: the evolution 
of the third, fourth and fifth closest points to the puncture, at 
r = {3, 4, 5}M/8192, are shown. We do not expect the closest 
two points to be reliable, due to finite-difference errors across 
the puncture, although they show similar behaviour, shown 
in the lower plot. 



the time a point takes to reach the inner horizon, 
parametrized by the point's isotropic coordinate r. The 
point dX r — M/2 is at the horizon at t = 0, and so 
"reaches" the inner horizon immediately. The time for 
points with r < M/2 to reach the inner horizon appears 
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FIG. 10: Time for grid points to pass the inner horizon, 
as a function of coordinate r. The innermost grid point is 
initially a,t R — 2031M, and reaches the inner horizon in 
T2M = 3.35M. 



to grow Hnearly as we move toward the puncture. Very 
close to the puncture, however, the time grows logarith- 
mically, and even the closest grid point, at r = M/8192 
and initially at i? = 2049M, reaches the inner horizon by 
about t = 3.3M. A similar result is shown in [M| . 



C. Late-time behaviour and approach to the 
stationary solution 

To follow the behaviour of R{r) at later times, we 
return to our standard convergence series simulations. 
Figure [TT] shows i? as a function of time for a grid 
point at r = 3M/32. We see that the point, having 
retreated quickly through the inner horizon, overshoots 
Rf) = 1.31Af (indicated by a dashed line in the figure), 
before returning and settling to a value just larger than 
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FIG. 12: Time development of K at the horizon, for simula- 
tions with initial lapse (a) a — 'ip~'^, and (b) a = 1. 



FIG. 11: The Schwarzschild coordinate i? of a point at r = 
3Af/32 as a function of time. The point overshoots Rq = 
1.31Af , indicated by a dashed horizontal line, but for t > 40M 
settles on a value slightly larger than Rq. 




0.0002 
0.0001 



-0.0002 
-0.000.? 




Low - Medium 
5.943-(Mcdium - High) 



Rq. We may expect that a point's Schwarzschild coordi- 
nate R can decrease but not increase: a point that falls 
into the black hole cannot rise back toward the surface, 
unless the lapse becomes negative and time progresses 
backward. We will return to this surprising behaviour 
when we represent the time development of the numeri- 
cal slices on a Penrose diagram. 

Having described the general behaviour of moving- 
puncture simulations of Schwarzschild wormhole punc- 
ture data, we would like to verify that the numerical so- 
lution converges to the analytic one. A natural quantity 
to look at would be R at the grid points closest to the ori- 
gin. In Figures 2 and 3 of [BT , we have seen that the met- 
ric quantities are approximately fourth-order convergent 
for at least 50M of evolution. However, the convergence 
is not so precise that it is retained in quantities derived 
from the evolution variables. In particular, R — '^'^XsJ'^^^ 
does not exhibit clean fourth-order convergence, and is 
not suitable for verifying convergence towards the ana- 
lytic 1-l-log geometry. 

As an alternative, we look at the value of Tr(i4r) on 
the horizon. Although locating the horizon once again 
requires an estimation of R, the accuracy is much better 
far from the origin, and a more systematic analysis of the 
convergence and accuracy of the solution is possible. 



Figure 12 shows the value of Tr(ii') at the horizon 
for simulations that begin with a = 1 and a = V""^- 
In both cases the value relaxes to the analytic result of 



Kh = 0.0668. Figure 13 shows the convergence of Kh for 
the convergence series described earlier. With a "precol- 
lapse" initial lapse of we see no sign of convergence 
at early times, but reasonably clean fourth-order conver- 
gence after about 15M of evolution. With an initial lapse 
of a = 1, we see fourth-order convergence after only a few 
M of evolution, although the convergence deteriorates at 
later times. 

Figure [14] shows the deviation of from the ana- 
lytic value on a logarithmic plot. The values are scaled 
assuming fourth-order convergence, and we see clearly 
that the disagreement between the numerical Kh and 
the value for the stationary 1-1- log slice converge to zero 
at fourth-order at late times. This provides strong evi- 



FIG. 13: Convergence of K at the horizon, for simulations 
with initial lapse (a) a = ip~^ ) (b) a = 1. The a — ip~^ 
simulations are not convergent at early times, while the a = 1 
simulations lose clean convergence after about 50M. 



dence that the numerical slice does indeed approach the 
analytic stationary slice with high accuracy. 



D. Penrose diagrams of numerical results 

We now represent the time development of the numer- 
ical slices on Penrose diagrams, using the technique de- 
scribed in Section [ill B[ Figure [T5[ shows such a diagram. 
The initial conditions are chosen such that the initial 
data are at T = and therefore correspond to the hori- 
zontal line between and i^j. During the evolution, the 
slices move upwards symmetrically in the diagram. We 
can clearly see that the points move quickly to the right 
as the slices move up, and very soon the region near i'j^ 
(the second asymptotically flat end) is extremely poorly 
resolved. In effect the numerical slices lose contact with 
the second asymptotically flat end and congregate near 
the cylinder at i? = 1.31M, which is shown by a dashed 
line. 



We can also see in Figure 15 (and as was also clear 
in Figure 111, that the slices flrst penetrate Rq, before 
retreating later to a location just outside Rq. We focus 
on this behaviour in Figure [T6[ 



Our initial reaction to Figure 11 might be that this is 
a numerical error: the Schwarzschild R associated with 
a point can decrease, but it cannot increase unless the 
lapse is negative and time progresses backwards. We have 
already seen that the lapse is everywhere non-negative, so 
this behaviour appears to be contradictory. However, the 
points can move to larger values of R and move forward in 



time with the aid of a non-zero shift. Figure 16 illustrates 
how this is possible. 

As a further illustration of this point, consider an ar- 
bitrary slice through the Schwarzschild solution. Choose 
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FIG. 14: Logarithm of the deviation of K at the horizon from 
the analytic value K = 0.066811. The deviations are scaled 
assuming fourth-order convergence, and the results indicate 
that K converges to the analytic value with fourth-order ac- 
curacy. 




FIG. 15: Penrose diagram produced from numerical data. 
We can clearly see the numerical slice retract from the sec- 
ond asymptotically flat end. The times shown are t/M = 
0,0.25,0.75,1.25,2.0,2.5,3.5,8. 



a — and = 1. The data points march along the 
slice, some with decreasing R, some with increasing R, 
depending on the slice we chose. Thus there is no con- 
nection between the change of R and the allowed time 
vectors. 

It should also be clear that if we were to run our simu- 
lation without a shift, then the slices would penetrate Rq, 
but would not be able to retreat to Rq at later times. This 
is illustrated in Figure [17] which was produced from a nu- 
merical simulation with = 0. This behaviour is in di- 
rect contrast to what happens in the case of true maximal 
slicing (where the maximal slicing condition is imposed 
throughout the evolution, and is not approached only 
asymptotically, as with the case of the maximal variant 
of the l-|-log condition, (21)), where the slices approach 
Rq = 3M/2, but cannot pass through it. Maximal slic- 
ing is an elliptic condition and thus generates "barriers" , 
while 1+log slicing is hyperbolic and so no barriers exist. 

Note once again that the slices are isometric across the 
throat, as shown in Figures 15 and 17 The shift only 




FIG. 16: A close-up of the region near the cylinder at 
Ro = 1.31M. Although the time has elapsed from ta ~ 4.5M 
to tt — ISAf (with non-negative lapse), R at the innermost 
gridpoint has increased, Ra < Rb- 




FIG. 17: Penrose diagram of an evolution identical to that 
used for Figure |15[ except that in this case the shift is zero 
throughout the evolution and the data points are joined in 
the plot. We see that the numerical slices no longer retract 
from the second asymptotically flat end, and now penetrate 
i?o, and stay there. The times shown are the same as in 
Figure [15] 



It is difficult to see what happens to the numerical 
points at late times in Figure [15] because all of the points 
bunch up in the upper right corner of the diagram. We 
can change this by choosing T < for the initial slice 
when constructing the diagram. Figure 18 shows the 



relabels points within the slices and grid points move 
accordingly, so the figures show identical slices covered 
by different numerical grids. 



slices at numerical times t — 30, 37.25, 40, 50M, with the 
initial time chosen as Schwarzschild time Tq — — 40M. 
We see clearly that the shces approach the cyUnder at 
R = 1.31M. In addition the figure shows the analytic 
solution evaluated at the same times, and we see that 
the numerical points lie perfectly on top of the analytic 
solution, and that since we have reached the stationary 
slice, the numerical and Schwarzschild times coincide. 



V. NUMERICAL EVOLUTION OF TRUMPET 
INITIAL DATA 

In the previous section we evolved wormholc puncture 
data with l-|-log slicing and the F-driver shift condition, 
and found that the numerical data quickly evolved from 
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FIG. 18: The numerical data at Schwarzschild times t = 
30, 37.25, 40, 50M, with the initial time chosen as Tq = -40M. 
The analytic solution, evaluated at the same times, is shown 
by a solid grey line. We see that the numerical points lie 
perfectly on the analytic solution. 
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FIG. 19: The error in -y^x (i.e., jxx — 1) after 50M of evo- 
lution of the time-independent trumpet puncture data. The 
errors for simulations at three resolutions are scaled assuming 
fourth-order convergence, and we indeed see that the errors 
converge to zero at fourth-order. 



a wormhole to a trumpet geometry. In this section we 
start with trumpet data. We first show explicitly that 
they are time-independent (up to numerical errors, which 
converge to zero with increasing resolution). We then 
demonstrate that, if we alternate between variants of 
l-}-log slicing during the evolution (in practice ( |21[ ) and 
( 22 1 ) , the numerical slice alternates accordingly between 
the respective stationary trumpet geometries. This pro- 
cess also allows us to illustrate how the coordinates can 
drift in these evolutions, while invariant quantities re- 
main unchanged. 
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A. Evolution of time-independent data 

If we evolve the stationary 1+log solution in isotropic 
coordinates, given in Section |IID 3| we expect the data 
to be time-independent. This is certainly the case when 
we look at the data by eye: the evolution variables do 
not appear to change at all. 

A more systematic test is shown in Figure |19[ where 
we show the error in j^x at t = 50M. The conformal 
spatial metric is flat in the stationary data, 'jij = 6ij, 
and so we can calculate the error by simply evaluating 
'^xx — 1- We see that some error has developed hy t — 
50M, and it has a peak at around x — 3M. The errors 
are scaled consistent with fourth-order convergence, and 
we indeed see fourth-order convergence up to around x = 
145M. Since the outer boundary is at a; = 192M, by 
t = 50M lower-order errors from the outer boundary will 
have propagated to around x = 142Af, and so we do not 
expect to see fourth-order convergence beyond this point. 

In Figure [20] we show the L2 norm of the error in 
along the x-axis as a function of time. The lines are 
once again scaled assuming fourth-order convergence to 
zero. We see reasonably clean fourth-order convergence, 
which appears to deteriorate slightly near the end of the 
simulation, although by this time lower-order errors from 



FIG. 20: The L2 norm of the error in ^y^x over the course of 
the evolution. The errors again converge to zero at fourth- 
order. 



the outer boundary will have contaminated the solution, 
as is clear in Figure [191 

These figures indicate that the data are indeed time 
independent, up to small numerical errors that converge 
to zero at the expected rate. Since the analytic value of 
"fxx is unity, we can easily calculate the percentage error 
from the figures: we see that a,t t — 50M the largest error 
in is 0.12%, for the lowest-resolution simulation. As 
such, we see that these data provide an excellent testing 
ground for the accuracy of a numerical code. This point 
needs to be emphasized. The moving-puncture approach 
as used here is currently the most popular method for 
simulating black-hole binaries. For such a code there is 
710 analytic black-hole solution that can be used to test 
the code, except for the 1-l-log and stationary maximal 
trumpet solutions presented here and in [33l [57] . These 
analytic solutions could be invaluable not only for test- 
ing a new code, but also in analysing and reducing the 
sources of error in current codes. 
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FIG. 21: The value of K at the horizon, for a simulation 
where the slicing condition alternates between standard 1+log 
and asymptotically maximal 1+log, Eqns. (36 1 and (21 1. The 



dashed horizontal lines indicate the respective analytic values 
of K at the horizon. 



FIG. 22: The coordinate location the horizon, for two simu- 
lations with an alternating slicing condition. The dashed line 
is for a simulation with the standard choice of ?7 = 2/M. The 
solid line is for a simulation with rj — 0. In this case, the 
horizon location almost returns to its original value when the 
solution returns to the 1+log stationary geometry. 



B. Alternating slices 

We start with 1+log trumpet data, as in the previous 
section. We evolve the data for t = 22.5M, and then 



switch the slicing evolution equation from (22 1 to (21 1, 
i.e., we change to the slicing condition that asymptotes 
to maximal slicing. After a further 50M of evolution, at 
t — 72.5M, we switch back to the original slicing condi- 
tion, which, assuming robustness of the method, should 
asymptote back to the stationary 1+log solution. 

Figure [21] shows the value of K at the horizon for this 
simulation. The value is K — 0.0668 on the horizon in the 
initial data, and remains at this value for the first 22. 5M 
of evolution. Then, when the slicing condition changes, 
K quickly evolves towards K = 0. Within about 30 Af 
we have « 0. At t = 72. 5M, the slicing condition 
is changed again, and within another 30M the slice has 
settled back to the stationary 1+log value. Simulations 
were performed with the same low, medium and high 



resolutions as the puncture data case in Section IV At 
t = 125M, when the simulations ended, the respective 
values of K on the horizon were 0.0679, 0.0670, 0.0669. 

These results illustrate the robustness of the moving- 
puncture method to locate the appropriate stationary 
1+log slice. Potentially more challenging tests have also 
been performed using excision initial data with the inte- 
rior filled in, and the method is again seen to be robust 

[SHllMllIQ]. 

In Figure [22] we show the coordinate distance r of 
the horizon from the origin. For the stationary 1+log 
data the horizon is at r = 0.8304M, and deviates by no 
more than 0.0014 % in the first 22M of evolution in the 
highest-resolution simulation. When the slicing condi- 
tion is switched to asymptotically maximal 1+log, and 
then back to standard 1+log, we see that the horizon 
does not return to its original position, at least not on 
the same time scale as the geometry returns to the sta- 
tionary 1+log geometry. (See the dashed line in the fig- 



ure.) This illustrates that, although the numerical slices 
quickly approach a stationary geometry, the coordinates 
may still drift. This effect is at least partially due to 
the damping parameter 77 in the F-driver shift condition. 
If we repeat the simulation with 77 = 0, we produce the 
solid line in Figure |22[ now the horizon location returns 
to its original location to a comparable accuracy that 
coordinate-invariant quantities return to the stationary 
1+log solution. 



VI. DISCUSSION 
A. Setting wormholes and trumpets in motion 

One of the directions for future work suggested by our 
results (and already proposed in Paper I) is the construc- 
tion of trumpet (as opposed to wormhole) initial data for 
black-hole binaries. 

Trumpet black-hole binary initial data would have the 
advantage that our numerical code would not need to 
evolve the fast transition from wormholes to trumpets. In 
addition, current binary puncture initial data start with 
zero speed across the numerical grid. When the simu- 
lation begins, the l+log/F-driver gauge conditions both 
transform the wormholes into trumpets, but also gener- 
ate an advection component to the shift vector, which 
moves the trumpets across the grid (Paper I). One could 
attempt to choose a "better" initial shift, so that the 
wormholes move from the outset, but if we begin with 
wormhole puncture data, there is no way to prevent this 
shift changing nontrivially as the wormholes evolve to 
trumpets. Although one could produce a "best" initial 
shift by some trial-and-error process, a more attractive 
proposal would be to start with trumpet data, and hope 
to find a procedure to choose an initial shift that imbues 
the boosted trumpets with their appropriate coordinate 
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speeds. We would then hope that, from the outset, almost 
all of the evolution of the data would represent the evolu- 
tion of physical quantities (the black holes' motions and 
the evolution of their spins, surface geometry, etc.) and 
mere gauge evolution would be minimised. This might 
help reduce the gauge component to the black-hole mo- 
tion seen in, for example, the first few hundred M of 
evolution in Figure 18 in |27j . Further, and more prac- 
tically, we would hope that such data would reduce the 
initial noise in wave quantities (for example. Figures 1, 
2, 4 and 6 in [27]). 

Before we become too optimistic, however, we should 
emphasize that the bulk of the initial noise in most nu- 
merical simulations probably comes from the burst of 
junk radiation present in the initial data. For example, 
the excision data evolved in [71, already have most of 
the properties we have just advertised: the gauge is not 
expected to evolve significantly in the early stages of the 
simulation, and the initial shift is precisely that which 
should cause the black holes to follow a quasi-circular in- 
spiral. As we would hope, the black-hole motion is indeed 
smooth from the outset. Nonetheless, noise still reduces 
the accuracy of some wave quantities at early times (see, 
for example, Figure 7 of [7J|), and the most likely culprit 
is the junk radiation, which has a similar magnitude to 
that in moving-puncture simulations. We expect that the 
ideal initial data for moving-puncture simulations would 
produce minimal junk radiation and be in trumpet form. 



B. Summary 

We have extended the analysis of the behaviour of the 
analytical and numerical slices in moving-puncture sim- 
ulations of the Schwarzschild spacetime that we began in 
Paper I. For a general form of the l-l-log slicing condition 
( 36 1 we have derived the stationary Schwarzschild trum- 



pet solution: the slice extends from spatial infinity to an 
infinitely long cylinder, or trumpet, with a throat at some 
finite radius Rq. When the parameter n in this condition 
is set to n = 2, we obtain the solution given in Paper I 
with Rq = 1.31AI . In the limit n — *■ cxd we recover the 
maximal trumpet solution f?31, '55], with Rq — 1.5M. 

For a given choice of the 1-1- log slicing condition, there 
is a unique regular stationary solution. In numerical sim- 
ulations that apply the moving-puncture technique to 
wormhole puncture Schwarzschild initial data, and use 
the F-driver shift condition, the numerical slice quickly 
evolves to the stationary slice. This cannot happen to 
the analytic slice: this must remain connected to the two 
asymptotically fiat ends in the wormhole data. It also 
cannot happen with numerical data if the shift is zero. 
However, the F-driver shift condition generates a shift 
that stretches the numerical slice such that all of the 
numerical points extremely quickly move onto the sta- 



tionary l-|-log slice; the non-stationary part of the slice 
no longer contains any grid points. The stretching of the 
slice is so extreme that no matter how many numerical 
points we place near the puncture (so long as there is 
no point on the puncture) that point will quickly move 
onto the stationary slice. Even a grid point initially at 
r = M/8192 and R « 2000A/ on the second copy of the 
exterior space, passes through R = 2M in less than 3.5M 
of evolution, and soon after settles near R — l.BlAf. 

An alternative to wormhole puncture data are trum- 
pet puncture data. We have transformed the stationary 
l-|-log solution to isotropic coordinates, and shown that 
these data are indeed time-independent when evolved, up 
to small numerical errors. In addition, we have shown 
that the numerical data can easily change from one sta- 
tionary geometry to another, if the l-|-log condition is 
changed during the evolution, indicating a certain ro- 
bustness of the method. Finally, we are able to see clearly 
that although the data approach a stationary slice, the 
numerical coordinates may still drift. 

The realization that the moving puncture method is 
really based on trumpet puncture data, and that this 
type of geometry naturally avoids most of the unphys- 
ical regions of spacetime in black hole evolutions, es- 
tablishes a new paradigm for the numerical evolution of 
black holes and suggests many directions of possible fu- 
ture research. One of the most promising is the con- 
struction of trumpet initial data for black-hole binaries, 
as discussed above. Additionally, one may use the sta- 
tionary data to make precise tests of numerical codes, 
for example to improve treatment of mesh-refinement 
boundaries and outer boundaries, to determine the res- 
olutions necessary to most accurately resolve black-hole 
spacetimes, and to explore different gauge choices. The 
stationary solution also provides an ideal background for 
mathematical studies of the stability properties of the 
BSSN/moving-puncture system and other evolution sys- 
tems, that could even be particularly tailored to evolve 
trumpet puncture data. 
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